1. Field of the Invention
The present invention relates to an imaging apparatus and method for displaying magnified images of sequentially estimated areas of current sources for estimating the position of an active area within a living body by using a biomedical magnetism measuring apparatus, including a highly sensitive magnetic field sensor, for example, a SQUID (Superconducting QUantum Interference Device) for measuring a magnetic field generated from a living body.
2. Description of the Related Art
With recent advances in superconductive device technology, highly sensitive biomedical magnetism measuring apparatus utilizing a SQUID have recently been employed in medical diagnostic apparatus. Such apparatus, which are also referred to as biomagnetometers, SQUID computed tomography (SQUID CT), magnetic source imaging (MSI), biomagnetic imaging (BMI), magnetoencephalograms (MEG), and magnetocardiograms (MEG) operate as follows. Electric sources within a living body simultaneously generate a low level magnetic field. Therefore, performing an inverse estimation (or inverse problem) of the active area within a living body by measuring the distribution of this magnetic field is expected to be useful for diagnosis of diseased regions within a living body. The terms inverse estimation and inverse problem refer to an algorithm in which a measured magnetic field is employed for estimating the location and/or distribution of the electric sources within the body. The position of a current source which acts as a magnetic field generating source is estimated to analyze heart disease and brain function disease from the measured magnetic field. For this purpose, a current dipole is used as a model of a biological source and is estimated within a homogeneous conductor using a computation model. A current dipole is a short segment of current which is used to illustrate a transient current flow in a small area. The inverse problem approach is employed for determining the current distribution for which the computed magnetic field becomes equal to the measured magnetic field. In accordance with the inverse problem approach, an algorithm is used to move the estimated value of the current dipole to the position for which the corresponding computed magnetic field approximates the measured magnetic field. The algorithm which is employed is based on the least square error solution (see Equation (5) below).
The above-described biomagnetic measuring apparatus should not be confused with magnetic resonance imaging (MRI) which detects only the configuration of a structure. Instead, the subject invention is directed to determining the functional state of an organ by detecting current paths in the body and particularly the brain (so-called neuromagnetism) and the heart (so-called cardiomagnetism). To indicate the magnitude of the magnetic fields created by such current flow in a living body, the neuromagnetic field is approximately 10.sup.-14 tesla, while the cardiomagnetic field is 10.sup.-12 tesla. The magnetic field is measured in order to determine the current amplitude and position of the equivalent current dipole.
Biomagnetometers are currently available in the art. For example, Biotechnology Incorporated (BTi) produces a neuromagnetometer utilizing a SQUID. Other biomagnetometers are manufactured by Siemens and CTF of Canada.
Referring to FIG. 18, in prior art biomedical magnetism measuring apparatus, a homogeneous semi-infinite conductor model of the torso, for example, is used for estimation of the current source in the heart (S1). Alternatively, a homogeneous or multilayer concentric conductor sphere can be used as a model for the head. The parameters of the current dipoles for the heart model are then estimated (S2) and the magnetic field B.sub.c based on these estimated dipoles is calculated (S3). The magnetic field of the heart in the living body is measured (S4) and the measured data B.sub.m is input to the computer (S5). A current dipole which minimizes an objective function equal to the squared difference between the measured magnetic field and the computed magnetic field, has been considered as the estimated position of the current source (also referred to above as the inverse problem or inverse estimation). Therefore, the measured value B.sub.m and the computed value B.sub.c are compared (S6) and if the squared error is minimized (S7) then the distribution of the current dipoles is displayed (S8). If the error is large (S7) then the parameters of the test dipoles are modified (S9) and the magnetic field B.sub.c is recalculated based on the newly set dipoles (S3).
The above-described approach for determining the amplitudes and positions of current dipoles is used in a number of different fields including the above-described biomedical magnetism apparatus and in other fields where the determination of current dipoles is desirable. However, certain problems are inherent in this approach. Local minima in the squared error between the measured and computed magnetic fields, can produce incorrect solutions to the dipole parameters. To overcome this problem, a long time has been required to obtain accurate current source localization because the computation does not converge within a finite amount of time to the correct solution (i.e., the global minimum of the objective function in the inverse problem).
As explained above, the least square error method for solving such a non-linear system requires repeated computation. As a method of avoiding such repeated computation, it has been suggested to fix the position of current dipoles on grid points, thereby making the problem linear. Such methods are described in, for example, Jeffs et al., "An Evaluation of Methods for Neuromagnetic Image Reconstruction", IEEE Transactions on Biomedical Engineering, Vol. BME-34, No. 9, September 1987, pp. 713-723; Smith et al., "Linear Estimation Theory Applied to the Reconstruction of a 3-D Vector Current Distribution", Applied Optics, Vol. 29, No. 5 (1990), pp. 658-667; and Sarvas, "Basic Mathematical and Electromagnetic Concepts of the Biomagnetic Inverse Problem", 1989, Phys. Med. Biol. 32, pp. 11-22, the contents of all which are incorporated by reference. A method for obtaining a least square solution by relating the measured magnetic field to the intensity of the current dipole by using simultaneous linear equations can be then formulated. To express the measured magnetic field obtained by pickup coils and the amplitude of the current dipole on each grid point with simultaneous linear equations, the distribution of the current dipoles is defined on a set of fixed grid points. The amplitudes in three directions of n current dipoles are defined as (q.sub.1x, q.sub.1y, q.sub.1z) . . . (q.sub.nx, q.sub.ny, q.sub.nz); the positions of the corresponding grid points are defined as (x.sub.1 ', y.sub.1 ', z.sub.1 ') (x.sub.n ', y.sub.n ', z.sub.n '); the amplitudes in three directions of the magnetic field measured by m pickup coils at m points are defined as (b.sub.1x, b.sub.1y, b.sub.1z) . . . (b.sub.mx, b.sub.my, b.sub.mz); the positions of the corresponding pickup coils are defined as (x.sub.z, y.sub.1, z.sub.1) . . . (x.sub.m, y.sub.m, z.sub.m), the current dipole vector Q={q.sub.1, q.sub.2, . . . q.sub.n }.sup.T and the measured magnetic field vector B.sub.m =(b.sub.1, b.sub.2, . . . b.sub.m).sup.T. Based on the Biot-Savart law, the simultaneous linear equation B=AQ can be solved, where the matrix of coefficients A is given by: ##EQU1##
Moreover, the elements of the equation B.sub.m =AQ can be expressed as follows: ##EQU2##
Here, from the Biot-Savart law, elements of each coefficient matrix can be obtained from the following expression: ##EQU3##
The equation B.sub.m =AQ is a linear equation which is determined by the current dipole positions and pickup coil positions. Therefore, the current dipole values Q can be obtained by solving this equation when m, the number of measurements, equals n, the number of unknowns. If the coefficient matrix A is a non-singular matrix, the inverse matrix A.sup.-1 exists and a current dipole distribution Q can be solved directly from EQU Q=A.sup.-1 B.sub.m ( 4)
However, if the coefficient matrix A is singular or if n&gt;m, the inverse matrix cannot be obtained and a unique solution does not exist. However, in this case, the product A.sup.T A of the coefficient matrix A and transposed matrix A.sup.T becomes a square matrix and A.sup.T A can be inverted when the column vectors of A are independent. In this case, the least squares solution, denoted Q.sub.c, which minimizes ##EQU4## which is the squared sum of the differences of measured values B.sub.m and computed values B.sub.c =AQ.sub.c, given by the normal equation EQU Q.sub.c =(A.sup.T A).sup.-1 A.sup.T B.sub.m ( 6)
can be obtained as described by Strang, "Linear Algebra and Its Applications", 1980, New York; ACADEMIC PRESS, INC., the contents of which is incorporated by reference. Under the least square error method of Equation (5), the parameter .beta. defines the test dipole position x', y' and z' and dipole strength q.sub.x, q.sub.y and q.sub.z. Equation (6) refers to the formulation which is one case of the least square error solution using singular value decomposition. Thus, the terms inverse estimation and inverse problem used above are based on the least square solution of Equation (5).
Moreover, when the matrix A is singular, where the column vectors of matrix A are not independent (i.e., rank (A)&lt;n), the inverse matrix of A.sup.T A does not exist and therefore a unique solution cannot be obtained. In this case, singular value decomposition can be utilized.
In accordance with singular value decomposition, a desired (m.times.n) matrix A can be decomposed to EQU A=U.LAMBDA.V.sup.T ( 7)
with the m.times.m orthogonal matrix U, m.times.n diagonal matrix .LAMBDA. and (n.times.n) orthogonal matrix V. .LAMBDA. is a diagonal matrix where the elements or singular values .lambda..sub.i (i=1, 2 . . . m) are the square roots of the eigenvalues of AA.sup.T and A.sup.T A which are arranged on the diagonal in descending order and U and V are eigenvectors of AA.sup.T and A.sup.T A, respectively, as described in Forsythe et al., "Computer Methods for Mathematical Computations, Prentice-Hall, New Jersey, (1978), the contents of which are incorporated by reference.
In this case, the least-square minimum-norm solution Q.sup.+ of Equation (4) above can be obtained from the following equation for the generalized inverse matrix EQU Q.sup.+ =V.LAMBDA..sup.+ U.sup.T B.sub.m =A.sup.+ B.sub.m ( 8)
Here, .LAMBDA..sup.+ is a diagonal matrix whose element is EQU .lambda..sup.+.sub.i =1/.lambda..sub.i ( 9)
when .lambda..sub.i is not equal to zero, .lambda..sup.+.sub.i =0 when .LAMBDA..sub.i =0. A.sup.+ is a pseudo inverse matrix, where the inverse matrix A.sup.-1 is extended to an arbitrary (m.times.n) matrix from an (n.times.n) square matrix.
A generalized inverse matrix method utilizing normal equations and a method utilizing singular value decomposition are effective for obtaining the density distribution of many current dipoles because a multi-dipole model is assumed. In this case, once the inverse matrix is obtained by using the normal equations or singular value decomposition, the current source density distribution Q can be obtained simply by multiplying the measured values B.sub.m either by the coefficient A.sup.+ as in Equation (8) or by (A.sup.T A).sup.-1 A.sup.T as in Equation (6). The distribution of current dipoles can then be obtained with higher speed than with the iterative method for solving the nonlinear least square system in which the position is estimated by moving current dipoles.
In the technique for expressing the relationship between magnetic field intensity and current dipole locations with simultaneous linear equations and obtaining a generalized inverse matrix from the normal equations, a method for dividing the higher current dipole intensity area based on the initial estimated value is known, and is described in Okada a al., "Current Density Imaging as a Method of Visualizing Neuronal Activity of the Brain", Society for Neuroscience Abstracts, 509.16:1241 (1990), the contents of which is incorporated by reference. However, in this method, current dipole resolution can be improved while leaving grid points existing in the peripheral area, but since the number of grid points of current dipole locations increases through division into subsections, the influence of the limited number of pickup coils limits the resolution. In addition, Okada's method assumes a current dipole plane (source plane) which is parallel to the measurement plane and does not describe a flexible grid point distribution method and a display method therefor.
As described above, the number of sensors m must be equal to or larger than the number of current dipoles n in order to accurately obtain the intensity distribution of current dipoles from the least square solution of such simultaneous linear equations. Therefore, a large number of pickup coils and SQUID magnetometers for measuring magnetic field intensity are required for estimating the position of current dipoles with an accuracy of several millimeters required for medical diagnosis. It has also been impossible to realize a real-time apparatus for displaying movement of dipole sources.